# Load necessary libraries
library(ggplot2)
library(Seurat)
library(tibble)
library(ggrepel)
library(dplyr)

#Load subsetted RPE
load("/nfs/turbo/umms-lprasov/Gabbi/RPEcombined_sub_seurat.Rdata")
RPEcombined_sub.seurat <- NormalizeData(RPEcombined_sub.seurat)

#Do differential gene expression and set logfc.threshold = 0 to get all of the genes
humRPE.markers <- FindMarkers(RPEcombined_sub.seurat, ident.1 = "RPE_Myrf Mutant", ident.2 = "RPE_Myrf Control", logfc.threshold = 0, verbose = FALSE)
P0RPE.markers <- FindMarkers(RPEcombined_sub.seurat, ident.1 = "RPE_P0-MUT", ident.2 = "RPE_P0_WT", logfc.threshold = 0, verbose = FALSE)

#Write files
write.csv(humRPE.markers, "humRPE.markers_SV5_021725.csv")
write.csv(P0RPE.markers, "P0RPE.markers_SV5_021725.csv")


#Rename columns
humRPE.markers$Geno <- "Hum_FC"
names(humRPE.markers)[names(humRPE.markers) == "avg_log2FC"] <- "HumFC"
P0RPE.markers$Geno <- "P0_FC"
names(P0RPE.markers)[names(P0RPE.markers) == "avg_log2FC"] <- "P0_FC"

#Remove any genes that are not present in both
RPE_merged_degs = merge (humRPE.markers, P0RPE.markers, by = 0, all = TRUE)
RPE_merged_degs_omit <- na.omit(RPE_merged_degs)
names(RPE_merged_degs_omit)[names(RPE_merged_degs_omit) == "avg_log2FC.x"] <- "HumFC"
names(RPE_merged_degs_omit)[names(RPE_merged_degs_omit) == "avg_log2FC.y"] <- "P0_FC"
RPE_merged_degs_omit1 <- RPE_merged_degs_omit[,-1]
rownames(RPE_merged_degs_omit1) <- RPE_merged_degs_omit[,1]

#Remove any genes with p-adjusted value greater than 0.1 in a gene for EITHER genotype...
RPE_merged_degs_omit2 <- RPE_merged_degs_omit1 %>% filter(RPE_merged_degs_omit1$p_val_adj.x <= 0.1 | RPE_merged_degs_omit1$p_val_adj.y <= 0.1)

RPE_merged_degs_omit1$P0_FC

RPE_merged_degs_omit11 <- RPE_merged_degs_omit1 %>% filter(abs(HumFC) > 0.25 & (abs(P0_FC) > 0.25))
RPE_merged_degs_omit_sigonly <- RPE_merged_degs_omit2 %>% filter(abs(HumFC) > 0.25 & (abs(P0_FC) > 0.25))

#Getting lists of genes in each category
RPE_mrgd_degs_sig_upup <- RPE_merged_degs_omit_sigonly %>% filter(HumFC > 0 & P0_FC > 0)
RPE_mrgd_degs_sig_downdown <- RPE_merged_degs_omit_sigonly %>% filter(HumFC < 0 & P0_FC < 0)
RPE_mrgd_degs_sig_Lup <- RPE_merged_degs_omit_sigonly %>% filter(HumFC > 0 & P0_FC < 0)
RPE_mrgd_degs_sig_Rdown <- RPE_merged_degs_omit_sigonly %>% filter(HumFC < 0 & P0_FC > 0)

RPE_mrgd_degs_sig_upup1 <- RPE_mrgd_degs_sig_upup %>% filter(p_val_adj.x <= 0.1 & p_val_adj.y <= 0.1)
RPE_mrgd_degs_sig_downdown1 <- RPE_mrgd_degs_sig_downdown %>% filter(p_val_adj.x <= 0.1 & p_val_adj.y <= 0.1)
RPE_mrgd_degs_sig_Lup1 <- RPE_mrgd_degs_sig_Lup %>% filter(p_val_adj.x <= 0.1 & p_val_adj.y <= 0.1)
RPE_mrgd_degs_sig_Rdown1 <- RPE_mrgd_degs_sig_Rdown %>% filter(p_val_adj.x <= 0.1 & p_val_adj.y <= 0.1)



RPE_merged_degs_omit11 <- RPE_merged_degs_omit11 %>% filter(RPE_merged_degs_omit11$p_val_adj.x <= 0.05 | RPE_merged_degs_omit11$p_val_adj.y <= 0.05)
                          

write.csv(RPE_merged_degs_omit2, "RPE_mergedDEGS_omit2.csv")


write.csv(RPE_merged_degs_omit2, "RPE_mergedDEGS_omit2.csv")

#Make the graph
RPE_merged_degs_omit1$gene <- rownames(RPE_merged_degs_omit1)
RPE_merged_degs_omit1$group <- ifelse(RPE_merged_degs_omit1$gene %in% c("Myrf"), "Humanized Gene",
                                      ifelse(RPE_merged_degs_omit1$gene %in% c("Wfikkn2", "Tcf4", "Id1", "Id3"), "TGFB/BMP", "Extracellular Matrix Marker"))

RPE_degs_plot2 <- ggplot(RPE_merged_degs_omit1, aes(x = P0_FC, y = HumFC)) +
  geom_smooth(method = "lm") +
  geom_point() +
  geom_text_repel(
    data = subset(RPE_merged_degs_omit1, gene %in% c("Myrf", "Wfikkn2", "Tcf4", "Id3", "Id1", "Upk3b", "Cdh3")),
    aes(label = gene, color = group),
    size = 5, fontface = "bold") +
  xlim(-3, 1.6) +
  ylim(-1.5, 1.5) +
  annotate("text", x = c(-2, -2), y = c(1.3, 1.1), label = c("p-value: < 2.2e-16", "Adjusted R-squared: 0.07154"), size = 7)+ theme_gray(base_size = 18) + guides(color = guide_legend(override.aes = list(shape = 16)))

#Make the graph WITHOUT THE CORRELATION LINE!!
RPE_merged_degs_omit_sigonly$gene <- rownames(RPE_merged_degs_omit_sigonly)
RPE_merged_degs_omit_sigonly$group <- ifelse(RPE_merged_degs_omit_sigonly$gene %in% c("Myrf"), "Humanized Gene",
                                      ifelse(RPE_merged_degs_omit_sigonly$gene %in% c("Wfikkn2", "Tcf4", "Id1", "Id3"), "TGFB/BMP", "Extracellular Matrix Marker"))

RPE_degs_plot2 <- ggplot(RPE_merged_degs_omit_sigonly, aes(x = P0_FC, y = HumFC)) +
  geom_point() +
  geom_hline(yintercept = 0)+
  geom_vline(xintercept = 0)+
  geom_label_repel(
    data = subset(RPE_merged_degs_omit_sigonly, gene %in% c("Myrf", "Wfikkn2", "Tcf4", "Id3", "Id1", "Upk3b", "Cdh3")),
    aes(label = gene, color = group),
    size = 5, fontface = "bold",
    box.padding = 1,   # Padding between box and point
    point.padding = 0.3,  # Padding between box and label
    segment.color = 'red',# Color of the segment line
    segment.size = 1
  ) +
  xlim(-9,5)+
  ylim(-10,3)


#Figure out how many genes are in each quadrant of the graph
#Positive concordant, UP in BOTH
sum((RPE_merged_degs_omit11['HumFC']>0) & (RPE_merged_degs_omit11['P0_FC']>0))
#Negative concordant, DOWN in BOTH
sum((RPE_merged_degs_omit11['HumFC']<0) & (RPE_merged_degs_omit11['P0_FC']<0))
#HumFC UP, P0_FC DOWN, Left upper quad
sum((RPE_merged_degs_omit11['HumFC']>0) & (RPE_merged_degs_omit11['P0_FC']<0))
#HumFC DOWN, P0_FC UP, Right bottom corner
sum((RPE_merged_degs_omit11['HumFC']<0) & (RPE_merged_degs_omit11['P0_FC']>0))

#Make the graph WITHOUT THE CORRELATION LINE!! Graph with only genes that have an adj. pval less than or equal to 0.1 for EITHER genotype
RPE_merged_degs_omit2$gene <- rownames(RPE_merged_degs_omit2)
RPE_merged_degs_omit2$group <- ifelse(RPE_merged_degs_omit2$gene %in% c("Myrf"), "Humanized Gene",
                                      ifelse(RPE_merged_degs_omit2$gene %in% c("Wfikkn2", "Tcf4", "Id1", "Id3"), "TGFB/BMP", "Extracellular Matrix Marker"))

RPE_degs_plot3 <- ggplot(RPE_merged_degs_omit2, aes(x = P0_FC, y = HumFC)) +
  geom_hline(yintercept = 0) + 
  geom_vline(xintercept = 0)+
  geom_point() +
  geom_text_repel(
    data = subset(RPE_merged_degs_omit2, gene %in% c("Myrf", "Wfikkn2", "Tcf4", "Id3", "Id1", "Upk3b", "Cdh3")),
    aes(label = gene, color = group),
    size = 5, fontface = "bold")+
  xlim(-10,5)+
  ylim(-10,5)

#Figure out how many genes are in each quadrant of the graph
#Positive concordant, UP in BOTH
sum((RPE_merged_degs_omit2['HumFC']>0) & (RPE_merged_degs_omit2['P0_FC']>0))
#Negative concordant, DOWN in BOTH
sum((RPE_merged_degs_omit2['HumFC']<0) & (RPE_merged_degs_omit2['P0_FC']<0))
#HumFC UP, P0_FC DOWN, Left upper quad
sum((RPE_merged_degs_omit2['HumFC']>0) & (RPE_merged_degs_omit2['P0_FC']<0))
#HumFC DOWN, P0_FC UP, Right bottom corner
sum((RPE_merged_degs_omit2['HumFC']<0) & (RPE_merged_degs_omit2['P0_FC']>0))

#Perform linear regression to get the Rsquared value and Pval
#http://127.0.0.1:10655/graphics/93393d72-2e04-48f1-88f0-d9c679489213.png
lmod = lm(HumFC~P0_FC, RPE_merged_degs_omit1)
summary(lmod)

#Linear regression results
#Call:
#  lm(formula = HumFC ~ P0_FC, data = RPE_merged_degs_omit2)

#Residuals:
#  Min       1Q   Median       3Q      Max 
#-1.26528 -0.09059 -0.00451  0.08887  1.70272 

#Coefficients:
#  Estimate Std. Error t value Pr(>|t|)    
#(Intercept) 0.020182   0.002269   8.893   <2e-16 ***
#  P0_FC       0.247540   0.010830  22.857   <2e-16 ***
  ---
#  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Residual standard error: 0.186 on 6766 degrees of freedom
#Multiple R-squared:  0.07168,	Adjusted R-squared:  0.07154 
#F-statistic: 522.4 on 1 and 6766 DF,  p-value: < 2.2e-16